// POWER: -snes_mf -pc_type lu -eps_nev 1 -eps_tol 1e-10 -nxyz 20 -o SECOND -eps_type power -jacobi 2
// ARNOLDI: -snes_mf -pc_type lu -eps_nev 1 -eps_tol 1e-10 -o SECOND -nxyz 20 -eps_type arpack -jacobi 2 -eps_ncv 10
// bin/neutronics_simple_eigen_test -d 2 -snes_mf -pc_type lu -eps_nev 1 -eps_tol 1e-4 -o SECOND -nxyz 20 -eps_type arpack -jacobi 2  -snes_rtol 1e-4 -snes_atol 1e-10 -snes_max_it 100 -newt -snes_mf_operator -eps_ncv 10
// ST: -eps_type power -st_type sinvert -eps_tol 1e-10 -st_pc_type none -st_matmode shell -eps_power_shift_type <constant or rayleigh> -st_shift 0.2344 

// POWER: bin/neutronics_simple_eigen_test -d 2 -snes_mf -pc_type lu -eps_nev 1 -eps_tol 1e-10 -nxyz 20 -o SECOND -eps_type power -jacobi 2
// ARNOLDI: bin/neutronics_simple_eigen_test -d 2 -snes_mf -pc_type lu -eps_nev 1 -eps_tol 1e-10 -o SECOND -nxyz 20 -eps_type arpack -jacobi 2 -eps_ncv 10

// POWER
// bin/neutronics_simple_eigen_test -d 1 -snes_mf -pc_type lu -eps_nev 1 -eps_tol 1e-10 -o SECOND -nx 200 -eps_type power -jacobi 2

// ARPACK
// bin/neutronics_simple_eigen_test -d 1 -snes_mf -pc_type lu -eps_nev 1 -eps_tol 1e-10 -o SECOND -nx 200 -eps_type arpack -jacobi 2 -eps_ncv 5

// ARPACK-NEWTON
// bin/neutronics_simple_eigen_test -d 1 -o SECOND -nx 200 -pc_type lu -eps_type arpack -eps_nev 1 -eps_tol 1e-4 -snes_rtol 1e-4 -snes_atol 1e-10 -snes_max_it 100 -newt -jacobi 2 -snes_mf_operator -eps_ncv 5


// IAEA
// bin/neutronics_simple_eigen_test -d 2 -snes_mf -pc_type lu -eps_nev 4 -eps_tol 1e-10 -eps_type arpack -jacobi 2 -eps_ncv 5

// bin/neutronics_simple_eigen_test -d 2 -pc_type lu -eps_type arpack -eps_nev 10 -eps_tol 1e-5 -snes_rtol 1e-6 -snes_atol 1e-12 -snes_max_it 100 -newt -jacobi 2 -snes_mf_operator -eps_ncv 20

// bin/neutronics_simple_eigen_test -d 2 -eps_type arpack -eps_nev 10 -eps_tol 5e-4 -newt -snes_atol 1e-11 -jacobi 0 -snes_converged_reason -kspprec_rtol 1e-4 -kspprec_pc_type sor -kspprec_pc_sor_omega 0.67 -kspprec_type preonly


// bin/neutronics_simple_eigen_test-devel -d 2 -st_pc_type icc -pc_type lu -st_pc_factor_shift_nonzero -pc_factor_shift_nonzero -eps_nev 5 -eps_type krylovschur -eps_tol 1e-4 -newt -no_ew -st_type sinvert -st_shift 1.1 -eps_power_shift_type rayleigh


#include "getpot.h"

#include "libmesh.h"
#include "mesh.h"
#include "mesh_refinement.h"
#include "mesh_generation.h"
#include "utility.h"

#include "neutronics_system.h"
#include "wrapper.h"
  //#include "karma_monitor.h"
#include "slepc_eigen_solver.h"

// Petsc Wrapper matrix and vector classes
#include "petsc_matrix.h"
#include "petsc_vector.h"
#include "vtk_io_vtk_indip.h"

#include "petscsnes.h"

#include "slepc_eigen_solver.h"

  //#include "petsc_snes_utils.h"
#include "problem.cc"

#include "slepceps.h"
  //#include "petscmat.h"

// Global objects
static char help[] = "";

// extern "C" PetscErrorCode Neutronics_eigen_jacobian (SNES , Vec x, Mat *jac, Mat *prec, MatStructure *msflag, void *) ;
// extern "C" PetscErrorCode Neutronics_eigen_residual (SNES snes, Vec x, Vec r, void *ctx) ;
PetscErrorCode orthonormalize_vector(Vec, Vec) ;

// int noperator = 0, prec_calls = 0, tot_operator_count = 0, tot_prec_count = 0 ;
Mat eigMatrix, Loss, Gain ;
double eigval;
bool matfree = false, is_prec_context = false ;
int prec_type ;	// for preconditioning type context
PetscBool prec_op ;

std::vector<Vec> evecs ;
std::vector<PetscReal> evals_r, evals_c ;
IP ip_ortho ;
std::vector<PetscBool> which_ortho ;
int total_ip_dots = 0 ;

// Begin the main program.
#undef __FUNCT__
#define __FUNCT__ "neutronics_simple_eigen_test"
int main (int argc, char** argv)
{
	EPS eps ;
	ST st ;
	const EPSType type ;
	
	//PetscMatrix<Real> *lossmat, *gainmat ;
	
	/// cpu-time related quantities
	double start, finish, run_time = 0;
	
	int refine = 0 ;

	PetscErrorCode 	ierr;
	bool iosoln = false, is_newton_ctx = false, test_jac = false ;
	
	bool is_unstructured=false, is_lobatto=false ;
	
	unsigned int dim = 1, psx = 5, psy = 5, psz = 5 ;	//	Dimension of problem, Number of Finite Elements
	std::string order = "FIRST" ;	// Order of FE basis function
	std::string family = "LAGRANGE";	// FE Family
	
	// Initialize KARMA, libMesh and any dependent libaries including PETSc
  LibMeshInit init (argc, argv);
  
  GetPot karma(argc, argv) ;
  
	{
		// Check for proper calling arguments.
		if (argc < 3)
		{
			std::cout << "\n\nHelp:\n" << help ;
			
			// invalid arguments ; print the help message and exit
			exit(3) ;
			
		}
		
		// Brief message to the user regarding the program name and command line arguments.
		// karma.print_options() ;
		{
			std::cout << "Running " << argv[0];

			for (int i=1; i<argc; i++)
				std::cout << " " << argv[i];
			std::cout << std::endl << std::endl;
		}
    
		// Read problem dimension from command line.  Use int instead of unsigned since the GetPot overload is ambiguous
		dim = karma("-d", 1) ;
		if (dim < 1 || dim > 3)
		{
			std::cout << "\nBy default, solving the 1-Dimensional problem\n" ;
			dim = 1 ;
		}
		std::cout << "Main: Dimension of problem = " << dim ;

		// Read number of elements from command line
		family = karma("-fam", "LAGRANGE") ;
		std::cout << "Main: FE Family for spatial discretization = " << family ;
		
		// Read FE order from command line
		order = karma("-o", "FIRST") ;
		std::cout << "Main: Order of spatial discretization = " << order ;
		
		// Read number of elements from command line
		psx = karma("-nxyz", 5) ;
		psy = psz = psx ;
		
		psx = karma("-nx", psx) ;
		psy = karma("-ny", psy) ;
		psz = karma("-nz", psz) ;
		
		// get the uniform refinement level needed
		refine = karma("-refine", 0) ;

		// get the type of treatment for the Loss operator
		prec_type = karma("-type", 0) ;
		
		prec_op = PETSC_FALSE ;
		if ( karma.search("-prec_op") )
			prec_op = PETSC_TRUE ;
				
		// flag for matrix-free computation
		matfree =  karma.search("-mfree") ;
				
		// flag for solving the eigen problem with newton method
		is_newton_ctx =  karma.search("-newt") ;
		
		// Decide if the solution is written to a gmv file
		iosoln = karma.search("-io") ;
	
		// Flag whether this is a restart run or not
		is_lobatto = karma.search("-lob") ;

		// Flag whether this is a restart run or not
		is_unstructured = karma.search("-unst") ;
		
		test_jac = karma.search("-test_jac") ;
	}
	
  
	// All data from command line have been read ; Now create and initialize objects for use
	{
    PerfLog perf_log ("Main");
		perf_log.push("create_physics");
    
    Problem1 problem ;
    
		// Creates a neutronics mesh
		Mesh mesh (problem.dimension) ;
		
#if (PROBLEMTYPE == 4)
		std::string mshFile = "input/2d_iaea.msh" ;
		mesh.read (mshFile) ;
#elif (PROBLEMTYPE == 5)
		std::string mshFile = "input/2d_bwr.msh" ;
		mesh.read (mshFile) ;
#elif (PROBLEMTYPE == 7)
		std::string mshFile = "input/2d_benchmark_anl.msh" ;
		mesh.read (mshFile) ;
#else		
      //mesh.discretize(psx, psy, psz, order, is_unstructured, true, is_lobatto) ;
    MeshTools::Generation::build_cube (mesh,
                                       psx,
                                       (dim>1) ? psy : 0,
                                       (dim>2) ? psz : 0,
                                       0, 1.,
                                       0, 1.,
                                       0, 1.,
                                       (dim==1)    ? EDGE2 : 
                                       ((dim == 2) ? QUAD4 : HEX8));
#endif
		PetscPrintf(PETSC_COMM_WORLD, "\nMesh read successfully.") ;
		
		if (refine > 0)
		{
			MeshRefinement mesh_refinement (mesh);

			mesh_refinement.uniformly_refine (refine);
		}
		
		// Now create our Neutronics-System physics object 
		NeutronicsSystem neutronics (mesh, problem, problem.Egroups(), order, family);
		
		// Let the Neutronics system know that we are solving the eigen-problem
		// Also indicate whether we are solving this in a matrix-free fashion
		// This call adds the loss and gain matrices to the system
		neutronics.initialize_eigen_context(matfree, NHEP) ;

		// Initialize the data structures for the equation system.
		neutronics.init() ;
		
		
#if (PROBLEMTYPE == 5)
		
		bool clean_mesh =  karma.search("-clean_mesh") ;
		if (clean_mesh)
		{
			// serial process
			MeshBase::const_element_iterator       el     = mesh.elements_begin();
			const MeshBase::const_element_iterator end_el = mesh.elements_end();
			
			std::vector<bool> is_void ;	is_void.resize(mesh.n_elem()) ;
			
			PetscPrintf(PETSC_COMM_WORLD, "\nTotal elements = %d", mesh.n_elem()) ;
			
			unsigned int ele_num = 0 ;
			for ( ; el != end_el; ++el)
			{
				// Store a pointer to the element we are currently working on. 
				const Elem* elem = *el;

				bool is_void_material = is_element_void_material(Egroups(), elem) ;
				// if (is_void_material)
					// PetscPrintf(PETSC_COMM_WORLD, "\nVoid: Elem = %d, Material = %d", elem->id(), elem->subdomain_id()) ;
			
				// const Parameters& params = elem->params() ;
				// PetscPrintf(PETSC_COMM_WORLD, "\nElem = %d, Material = %d, Surface = %d", elem->id(), 
						// params->get<unsigned int>("physical"), params->get<unsigned int>("elementary")) ;
			
				is_void[elem->id()] = is_void_material ;
				
				ele_num++ ;
			}
			
			PetscPrintf(PETSC_COMM_WORLD, "\nNow deleting void elements") ;
			for (int i=0; i<is_void.size();  i++)
			{
				// Store a pointer to the element we are currently working on. 
				if (!is_void[i])
					continue ;
				
				mesh.delete_elem(mesh.elem(i)) ;
			}
			
			// since we deleted the elements, lets renumber the ids for elems and nodes
			mesh.renumber_nodes_and_elements () ;
			
			PetscPrintf(PETSC_COMM_WORLD, "\nWriting new mesh to file: 2d_bwr_modified.msh") ;
			GmshIO(mesh).write ("2d_bwr_modified.msh");
			std::cin.get() ;
		}
#endif
		
		perf_log.pop("create_physics");

		neutronics.print_info () ;
		
		// Now let's reset the snes context function pointers to the left and right multiply functions 
		perf_log.push("gain_loss_create");

		// Create our linear system; Allocate the loss and production matrices here
		{
			ierr = neutronics.assemble_gain_matrix () ; CHKERRQ(ierr);

			ierr = neutronics.assemble_loss_matrix () ; CHKERRQ(ierr);
		}

		Gain = neutronics.get_gain_matrix()->mat() ;
		Loss = neutronics.get_loss_matrix()->mat() ;
		
		perf_log.pop("gain_loss_create");
		
		PetscVector<Real>& pT = neutronics.get_residual () ;
		PetscVector<Real>& pS = neutronics.get_solution () ;
		
		int numdofs = pT.size() ;
		int locdofs = pT.local_size () ;
		
		if (test_jac)
		{
		
			numdofs += 1 ;
			locdofs = numdofs ;
			
			PetscVector<Real> temp1(numdofs, locdofs) ;
			PetscVector<Real> temp2(numdofs, locdofs) ;
			PetscVector<Real> temp3(numdofs, locdofs) ;
			PetscVector<Real> temp4(numdofs, locdofs) ;
			PetscVector<Real> temp5(numdofs, locdofs) ;
			PetscVector<Real> temp6(numdofs, locdofs) ;
			PetscVector<Real> temp7(numdofs, locdofs) ;
			
			Wrapper::initialize_wrapper(pS) ;
		
			double epsil = 1e-7 ;
			ierr = PetscOptionsGetReal(PETSC_NULL, "-test_jac_eps",&epsil,PETSC_NULL);
	
			std::cout << "about to start calculations: checking jacobian definition." ;
			
			// set a random vector and a perturbation vector
			ierr = VecSetRandom(temp1.vec(), PETSC_NULL) ;CHKERRQ(ierr);	// x
			ierr = VecSetRandom(temp3.vec(), PETSC_NULL) ;CHKERRQ(ierr);	// v
			
			std::cout << "solution and random vectors initialized." ;
			
			// set the solution to the random vector
			ierr = Wrapper::petsc_copy_vector(temp1.vec(), pS.vec(), false) ;CHKERRQ(ierr);
			neutronics.Keff() = temp1(numdofs-1) ;
			
			std::cout << "copied solution and updated keff." ;
			
			// 
			ierr = Wrapper::simple_eigen_residual(&neutronics, temp1.vec(), temp2.vec()) ;CHKERRQ(ierr);
			ierr = VecCopy(temp3.vec(), temp5.vec()) ;CHKERRQ(ierr);
			ierr = VecAYPX(temp3.vec(), epsil, temp1.vec()) ;	CHKERRQ(ierr);	// x+epsil*v
			ierr = Wrapper::simple_eigen_residual(&neutronics, temp3.vec(), temp4.vec()) ;CHKERRQ(ierr);
			
			std::cout << "solution and random vectors initialized." ;
			
			// temp4 = f (x+epsil*v) - f(x)
			ierr = VecAXPY(temp4.vec(), -1.0, temp2.vec()) ;
			// temp4 = (f (x+epsil*v) - f(x)) / epsil = J*v
			ierr = VecScale(temp4.vec(), 1.0/epsil) ;
			// Now temp4 holds the J*v vector
			
			Mat Jac ;
			ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)&neutronics, &Jac);CHKERRQ(ierr);
			ierr = MatSetFromOptions(Jac) ;CHKERRQ(ierr);
			
			ierr = MatShellSetOperation(Jac, MATOP_MULT,(void(*)())Wrapper::multiply_true_jacobian_matrix);CHKERRQ(ierr);
			ierr = MatShellSetOperation(Jac, MATOP_DUPLICATE,(void(*)())Wrapper::shell_matrix_duplicate);CHKERRQ(ierr);
		
			ierr = MatMult(Jac, temp5.vec(), temp6.vec()) ;CHKERRQ(ierr);
			
			ierr = VecCopy(temp6.vec(), temp7.vec()) ;
			// error in J*v from F.D perturbation vs Analytical jacobian
			// truly, this should be of O(epsil). If not, there is a problem !
			ierr = VecAXPY(temp7.vec(), -1.0, temp4.vec()) ;
			
			PetscScalar *fd_arr, *anal_arr, *err_arr ;
			ierr = VecGetArray(temp4.vec(), &fd_arr) ;
			ierr = VecGetArray(temp6.vec(), &anal_arr) ;
			ierr = VecGetArray(temp7.vec(), &err_arr) ;
			
			for (unsigned int i=0; i < temp4.size(); i++)
			{
				PetscPrintf(PETSC_COMM_WORLD, "\nF.D=%G\tANAL=%G\tERROR=%G", fd_arr[i], anal_arr[i], err_arr[i]) ;
			}
			
			Real errnrm = temp7.l2_norm() ;
			PetscPrintf(PETSC_COMM_WORLD, "\nNorm of error = %G", errnrm) ;
			
			ierr = VecRestoreArray(temp4.vec(), &fd_arr) ;
			ierr = VecRestoreArray(temp6.vec(), &anal_arr) ;
			ierr = VecRestoreArray(temp7.vec(), &err_arr) ;
			
			Wrapper::clear_wrapper() ;
			
			exit(3) ;
		}
		
		std::cout << "done." ;

		// Create eigensolver context
		ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr);
		
		ierr = EPSGetST(eps, &st) ;CHKERRQ(ierr);
		
		// Create the shell matrix
		{
			ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)&neutronics, &eigMatrix);CHKERRQ(ierr);
			ierr = MatSetFromOptions(eigMatrix) ;CHKERRQ(ierr);
			
			// ierr = MatShellSetOperation(eigMatrix, MATOP_MULT,(void(*)())Eig_Mult_Vector);CHKERRQ(ierr);
			ierr = MatShellSetOperation(eigMatrix, MATOP_MULT, (void(*)())Wrapper::simple_eigen_operator);CHKERRQ(ierr);
			ierr = MatShellSetOperation(eigMatrix, MATOP_SHIFT, (void(*)())Wrapper::shell_matrix_apply_shift);CHKERRQ(ierr);
			ierr = MatShellSetOperation(eigMatrix, MATOP_COPY,(void(*)())Wrapper::shell_matrix_copy);CHKERRQ(ierr);
			ierr = MatShellSetOperation(eigMatrix, MATOP_DUPLICATE, (void(*)())Wrapper::shell_matrix_duplicate);CHKERRQ(ierr);
			ierr = MatShellSetOperation(eigMatrix, MATOP_ZERO_ENTRIES,(void(*)())Wrapper::shell_matrix_zero_entries);CHKERRQ(ierr);
			
		}

		Wrapper::initialize_wrapper(pS.vec()) ;
		
		// set properties for ST in EPS
		{
			ierr = EPSSetOperators(eps, eigMatrix, PETSC_NULL) ;CHKERRQ(ierr);
			ierr = EPSSetProblemType(eps, EPS_NHEP) ;CHKERRQ(ierr);
			// ierr = EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE) ;CHKERRQ(ierr);
			
			// Set solver parameters at runtime
			ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
			
			// first check if shifted-inverse mode is ON. 
			const STType st_type ;
			ierr = STGetType(st, &st_type) ;CHKERRQ(ierr);

			if (!strcmp(st_type, "sinvert"))
			{
				// // we are using the shifted-inverse operator
				
				// set mat-mode to shell if other than constant shift is used
				// for some reason, inplace does not go beyond 1 iteration with Rayleigh and Wilkinson !!
				// EPSPowerShiftType pow_shift_type ;
				// ierr = EPSPowerGetShiftType(eps, &pow_shift_type) ;CHKERRQ(ierr);
				// if (pow_shift_type != 0)
				// {
					// std::cout << "Setting ST-MatMode to Shell." ;
					// ierr = STSetMatMode(st, STMATMODE_SHELL ) ;CHKERRQ(ierr);
				// }
			
				// perform few power iterations to get the shift
				
				double user_shift = karma("-st_shift", 0.0) ;	// by default - no preconditioner
				
				std::string pc_option ;     // preconditioner type - if present
				pc_option = karma("-st_pc_type", "none") ;	// by default - no preconditioner

				bool shell_pc = false ;
				if (!pc_option.compare("shell"))
					shell_pc = true ;
			
				// // now set mat mode to shell or else shifted-inverse based iterations do not converge
				KSP st_ksp ;
				ierr = STGetKSP(st, &st_ksp ) ;CHKERRQ(ierr);
				PC st_pc ;
				ierr = KSPGetPC(st_ksp, &st_pc) ;CHKERRQ(ierr);
				
				if (shell_pc)
				{
					// set our shell PC for ST matrix
					ierr = PCSetType(st_pc, PCSHELL) ;CHKERRQ(ierr);
					PCShellSetContext(st_pc, (void*)&neutronics) ;
					PCShellSetApply(st_pc, (PetscErrorCode (*)(PC,Vec,Vec))Wrapper::simple_eigen_preconditioner_solve_operator) ;
				}
				else
				{
          std::cout << "Specified -st_pc_type " << pc_option ;
					std::cout << "-st_pc_type shell not speicifed. Hence NO-PC will be used for ST solves.\n" ;
					
					// // set the st_pc_type = none since we do not have any information about shifted operators
					// // and hence cannot provide preconditioners; and of course, this is a shell matrix which
					// // makes the whole point of preconditioning even harder ! bah ...
					ierr = PCSetType(st_pc, PCNONE) ;CHKERRQ(ierr);
				}
				
				// by default - do not use this to jump-start sinvert method
				bool init_power_shift = karma.search("-init_power_shift") ;
				
				if (init_power_shift)
				{
					//////////////////////////////////
					double init_shift = karma("-init_shift", 2.5) ;	// by default - no preconditioner
					double init_tol = karma("-init_tol", 1e-2) ;	// by default - no preconditioner
					int init_maxit = karma("-init_maxit", 20) ;	// by default - no preconditioner
					double orig_tol ;
					int orig_maxit ;
					EPSGetTolerances(eps,&orig_tol,&orig_maxit);
					STSetShift(st, init_shift) ;	// no shift - pure power
					int nev_req, ncv_req ;
#if SLEPC_VERSION_LESS_THAN(3,0,0)
					ierr = EPSGetDimensions(eps, &nev_req, &ncv_req) ;CHKERRQ(ierr);
					ierr = EPSSetDimensions(eps, 1, 1) ;CHKERRQ(ierr);
#else
					ierr = EPSGetDimensions(eps, &nev_req, &ncv_req, 0) ;CHKERRQ(ierr);
					ierr = EPSSetDimensions(eps, 1, 1, PETSC_DECIDE) ;CHKERRQ(ierr);
#endif
					EPSSetTolerances(eps, init_tol, init_maxit) ;	// maximum 20 iterations to get first 2 digits
					
					// Solve the eigensystem
					ierr = EPSSolve(eps);CHKERRQ(ierr);
				
					int init_converged = 0 ;
					EPSGetConverged(eps, &init_converged) ;

					double eig_guess = 0 ;
					if (init_converged > 0)
					{
						ierr = EPSGetEigenpair(eps, 0, &eig_guess, PETSC_NULL, pS.vec(), PETSC_NULL);CHKERRQ(ierr);
					}
					
					eig_guess += user_shift ;
					//ierr = EPSSetInitialVector(eps, pS.vec()) ;CHKERRQ(ierr);
					ierr = EPSSetType(eps, EPSPOWER) ;CHKERRQ(ierr);
					ierr = EPSSetTolerances(eps, orig_tol, orig_maxit) ;CHKERRQ(ierr);
#if SLEPC_VERSION_LESS_THAN(3,0,0)
					ierr = EPSSetDimensions(eps, nev_req, ncv_req) ;CHKERRQ(ierr);
#else
					ierr = EPSSetDimensions(eps, nev_req, ncv_req, PETSC_DECIDE) ;CHKERRQ(ierr);
#endif
					ierr = STSetShift(st, eig_guess) ;	CHKERRQ(ierr);
					
					//////////////////////////////////////
				}
				
			}
		}
		
		int maxit ;
		double tol ;
		
		// Optional: Get some information from the solver and display it
		EPSGetType(eps,&type);
		EPSGetTolerances(eps,&tol,&maxit);
		int nev_req, ncv_req ;
#if SLEPC_VERSION_LESS_THAN(3,0,0)
		ierr = EPSGetDimensions(eps, &nev_req, &ncv_req) ;CHKERRQ(ierr);
#else
		ierr = EPSGetDimensions(eps, &nev_req, &ncv_req, 0) ;CHKERRQ(ierr);
#endif

		perf_log.push("compute_eigen_modes");
		start = MPI_Wtime();
		int nconv ;
		
		if(!is_newton_ctx)
		{

			// Optional: Get some information from the solver and display it
			PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
			PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev_req);
			PetscPrintf(PETSC_COMM_WORLD," Maximum dimension of the subspace: %d\n",ncv_req);
			PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);
			
			// Solve the eigensystem
			ierr = EPSSolve(eps);CHKERRQ(ierr);
			
			// {
				// SlepcEigenSolver<Real> eig_solver ;
				// eig_solver.set_eigensolver_type(KRYLOVSCHUR) ;
				// eig_solver.set_eigenproblem_type(GNHEP) ;
				// eig_solver.set_position_of_spectrum(LARGEST_MAGNITUDE) ;

				// std::pair<unsigned int, unsigned int> eig_result = eig_solver.solve_generalized(*gainmat,
								   // *lossmat,
								   // 1,	// we want only fundamental mode
								   // 40,	// set default subspace size
								   // tol,	// can make this user specified later
								   // maxit);	// for Krylov-Schur, this should be plenty.
				// ierr = EPSDestroy(eps) ;
				// eps = eig_solver.eps() ;
			// }
			
			int its ;
			ierr = EPSGetIterationNumber(eps, &its);CHKERRQ(ierr);
			
			int ops, dots, lits ;
			ierr = EPSGetOperationCounters(eps, &ops, &dots, &lits); CHKERRQ(ierr);
			PetscPrintf(PETSC_COMM_WORLD," Operator applications: %d\n",ops);
			PetscPrintf(PETSC_COMM_WORLD," Number of inner products: %d\n",dots);
			PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations: %d\n",lits);
			PetscPrintf(PETSC_COMM_WORLD,"\nNumber of iterations of the method: %d\n",its);
			
			// Display solution and clean up
			
			// Get number of converged approximate eigenpairs
			EPSGetConverged(eps,&nconv);
			PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %d\n\n",nconv);
			
			if (nconv>0)
			{
				PetscScalar kr, ki ;
	
				evecs.resize(nconv) ;
				evals_r.resize(nconv) ;
				evals_c.resize(nconv) ;
				
				for(int i=0; i<nconv; i++ )
				{
					ierr = VecDuplicate(pS.vec(), &evecs[i]) ;CHKERRQ(ierr);
					
					// Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
					// ki (imaginary part)
					ierr = EPSGetEigenpair(eps,i,&kr,&ki,evecs[i],PETSC_NULL);CHKERRQ(ierr);
					evals_r[i] = kr ;	evals_c[i] = ki ;
					
					// ierr = EPSGetRightVector(eps,i,evecs[i],PETSC_NULL);CHKERRQ(ierr);
				}
			}
		}
		else
		{
			numdofs += 1 ;
      locdofs = numdofs ;
			
			PetscVector<Real> eigres(numdofs, locdofs) ;
			PetscVector<Real> eigsol(numdofs, locdofs) ;
	
			perf_log.push("get_initial_guess");
			PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
			PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev_req);
			PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);
			
			// Solve the eigensystem to coarse tolerance using EPS object
			ierr = EPSSolve(eps);CHKERRQ(ierr);
			
			int its ;
			ierr = EPSGetIterationNumber(eps, &its);CHKERRQ(ierr);
			int ops, dots, lits ;
			ierr = EPSGetOperationCounters(eps, &ops, &dots, &lits); CHKERRQ(ierr);
			PetscPrintf(PETSC_COMM_WORLD," Operator applications: %d\n", ops);
			PetscPrintf(PETSC_COMM_WORLD," Number of inner products: %d\n", dots);
			PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations: %d\n", lits);
			PetscPrintf(PETSC_COMM_WORLD,"\nNumber of %s iterations: %d\n", type, its);
			
			perf_log.pop("get_initial_guess");
			
			Mat Jac, Prec ;
			{
				ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)&neutronics, &Jac);CHKERRQ(ierr);
				ierr = MatSetFromOptions(Jac) ;CHKERRQ(ierr);
				
				ierr = MatShellSetOperation(Jac, MATOP_MULT,(void(*)())Wrapper::multiply_true_jacobian_matrix);CHKERRQ(ierr);
				ierr = MatShellSetOperation(Jac, MATOP_SOLVE,(void(*)())Wrapper::preconditioner_solve);CHKERRQ(ierr);
				ierr = MatShellSetOperation(Jac, MATOP_DUPLICATE,(void(*)())Wrapper::shell_matrix_duplicate);CHKERRQ(ierr);
			
				ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)&neutronics, &Prec);CHKERRQ(ierr);
				ierr = MatSetFromOptions(Prec) ;CHKERRQ(ierr);
				
				ierr = MatShellSetOperation(Prec, MATOP_MULT,(void(*)())Wrapper::multiply_true_jacobian_matrix);CHKERRQ(ierr);
				ierr = MatShellSetOperation(Prec, MATOP_SOLVE,(void(*)())Wrapper::preconditioner_solve);CHKERRQ(ierr);
				// ierr = MatShellSetOperation(Prec, MATOP_MULT,(void(*)())Wrapper::preconditioner_multiply);CHKERRQ(ierr);
				ierr = MatShellSetOperation(Prec, MATOP_DUPLICATE,(void(*)())Wrapper::shell_matrix_duplicate);CHKERRQ(ierr);
			
			}
			
			SNES snes ;
			ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);
			
			ierr = SNESMonitorSet(snes, snes_monitor, &neutronics, 0);CHKERRQ(ierr);

			PetscBool ksp_ew_flag ;
			ierr = PetscOptionsHasName(PETSC_NULL,"-no_ew", &ksp_ew_flag);CHKERRQ(ierr);
			if (!ksp_ew_flag)
			{
				ierr = SNESKSPSetUseEW(snes, PETSC_TRUE) ;CHKERRQ(ierr);
				ierr = SNESKSPSetParametersEW(snes, 3, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT,
									PETSC_DEFAULT, PETSC_DEFAULT) ;
			}
			
			// In case of -sinvert mode, let us reset the shift and start clean
			// But for safety, let just do this all the time !
			neutronics.eigen_shift() = 0.0 ;
			
			ierr = SNESSetFunction(snes, eigres.vec(), eigen_residual, &neutronics);CHKERRQ(ierr);
			
			// bool use_preconditioner = karma.search("-init_power_shift") ;
			ierr = SNESSetJacobian(snes, Jac, Prec, eigen_jacobian, &neutronics); CHKERRQ(ierr);
			ierr = SNESSetTolerances(snes, 1e-10, 1e-5, PETSC_DEFAULT, 1000, 2000) ; CHKERRQ(ierr);
			
			KSP _ksp; PC _pc ;
			ierr = SNESGetKSP(snes, &_ksp) ;CHKERRQ(ierr);
			ierr = KSPGetPC(_ksp, &_pc) ;CHKERRQ(ierr);
			ierr = KSPSetType(_ksp, KSPFGMRES);CHKERRQ(ierr);
			
			bool init_knoll = karma.search("-init_knoll") ;
			if (init_knoll)
			{
				ierr = KSPSetInitialGuessKnoll(_ksp, PETSC_TRUE) ;CHKERRQ(ierr);
			}
			
			ierr = KSPSetPCSide(_ksp, PC_RIGHT) ;CHKERRQ(ierr);
			
			ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
			
			// always use "basic" newton iteration
			bool use_linesearch = karma.search("-cubic") ;
			if (!use_linesearch)
			{
				ierr = SNESLineSearchSet(snes, SNESLineSearchNo, PETSC_NULL) ; CHKERRQ(ierr);
			}
			
			PetscBool preconditioner_flag ;
			ierr = PetscOptionsHasName(PETSC_NULL,"-precond", &preconditioner_flag);CHKERRQ(ierr);
			if (preconditioner_flag)
			{	
				ierr = PCSetType(_pc, PCSHELL) ;
				PCShellSetContext(_pc, (void*)&neutronics) ;
				PCShellSetApply(_pc, (PetscErrorCode (*)(PC,Vec,Vec))Wrapper::preconditioner_solve) ;
			}
			else
			{
				ierr = PCSetType(_pc, PCNONE) ; CHKERRQ(ierr);
			}
			
			perf_log.push("hybrid_newton_solve");
			
			// Create an Orthogonalization context
			ierr = IPCreate(PETSC_COMM_WORLD, &ip_ortho) ;CHKERRQ(ierr);
			ierr = IPSetOptionsPrefix(ip_ortho, "newt_") ;
			ierr = IPSetOrthogonalization(ip_ortho, IP_ORTHOG_CGS, IP_ORTHOG_REFINE_IFNEEDED, 1e-3) ;CHKERRQ(ierr);
			
			ierr = IPSetFromOptions(ip_ortho) ;CHKERRQ(ierr);
			
        //Wrapper::set_orthonormalize_reference((PetscErrorCode(*)(Vec, Vec))orthonormalize_vector) ; 
			
			// std::vector<PetscBool> which_ortho ;
			// PetscScalar *ortho_h, *swork ;
			// PetscReal ortho_norm ;
			// PetscBool lindep ;
			// int ndots ;
			
			double kr, ki ;
			nconv = 0 ;
			// tot_operator_count = noperator ;	tot_prec_count = 0 ;
			for (int iev = 0 ; iev < nev_req; iev++)
			{
				ierr = EPSGetEigenpair(eps,iev,&kr,&ki,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
				ierr = EPSGetEigenvector(eps,iev,pS.vec(),PETSC_NULL);CHKERRQ(ierr);
				
				// // let us orthogonalize the input vector against all the previously converged vectors
				// ierr = IPOrthogonalize(ip_ortho, nconv, &which_ortho[0], &evecs[0], pS.vec(), 
								// ortho_h, &ortho_norm, &lindep, pT.vec()) ;CHKERRQ(ierr);
	
				// ierr = IPView(ip_ortho, PETSC_VIEWER_STDOUT_WORLD) ;CHKERRQ(ierr);
				// ierr = IPGetOperationCounters(ip_ortho, &ndots) ;CHKERRQ(ierr);
				// // if (lindep == PETSC_FALSE)
					// // PetscPrintf(PETSC_COMM_WORLD,"IPOrthogonalize refinement did not improve the quality of orthogonalization.\n");
				// PetscPrintf(PETSC_COMM_WORLD,"Number of inner-products during orthogonalization: %d",ndots);

				Real l2nrm = pS.l2_norm() ;
				ierr = VecNormalize(pS.vec(), &l2nrm) ;
								
				ierr = Wrapper::petsc_copy_vector(pS.vec(), eigsol.vec(), true) ;CHKERRQ(ierr);
				ierr = VecSetValue(eigsol.vec(), numdofs-1, kr, INSERT_VALUES) ;CHKERRQ(ierr);

				// // use P-1*res = -del*soln
				// ierr = neutronics_eigen_residual(snes, eigsol.vec(), eigres.vec(), &neutronics);CHKERRQ(ierr);
				// ierr = VecScale(eigres.vec(), -1.0) ;
				// Wrapper::preconditioner_solve(&neutronics, eigres.vec(), eigsol.vec()) ;
				// ierr = VecCopy(pT.vec(), pS.vec()) ;
								
				// set the keff value in neutronics system
				neutronics.Keff() = kr ;
				
				//ierr = VecCopy(eps->V[iev], ctx->x) ;CHKERRQ(ierr);
				//ierr = EPSGetRightVector(eps,0,ctx->x,PETSC_NULL);CHKERRQ(ierr);
				//ierr = EPSGetLeftVector(eps,0,temp,PETSC_NULL);CHKERRQ(ierr);
				
				// VecView(eigsol.vec(), 	PETSC_VIEWER_STDOUT_WORLD ) ;
				// std::cin.get() ;
				
				// calculate and set the initial S.S value ; could do a nonlinear solve to find this if exact function is unavailable in non-test mode
				//ierr = neutronics.form_initial_guess(*ctx->px);CHKERRQ(ierr);
			
				// prec_calls = 0 ;
				// noperator = 0 ;	// Start with number of Operator applications for EPS Solve
				// Solve the initial S.S problem
				ierr = SNESSolve(snes, PETSC_NULL, eigsol.vec()); CHKERRQ(ierr);
				
				SNESConvergedReason newton_convergence_reason;
				ierr = SNESGetConvergedReason(snes, &newton_convergence_reason) ; CHKERRQ(ierr);

				if(newton_convergence_reason >= 0)
				{
					evecs.resize(nconv+1) ;	evals_r.resize(nconv+1) ;	evals_c.resize(nconv+1) ;
					
					ierr = VecDuplicate(pT.vec(), &evecs[nconv]) ;CHKERRQ(ierr);
					ierr = Wrapper::petsc_copy_vector(eigsol.vec(), evecs[nconv], false) ;CHKERRQ(ierr);
					
					evals_r[nconv] = eigsol(numdofs-1) ;
					
					int lits, nits ;
					ierr = SNESGetIterationNumber(snes, &nits) ;CHKERRQ(ierr);
					ierr = SNESGetLinearSolveIterations(snes, &lits) ;CHKERRQ(ierr);
					
					// Get some statistics for the convergence of the nonlinear system
					PetscPrintf(PETSC_COMM_WORLD,"\nNumber of iterations of the method: %d",nits);
					PetscPrintf(PETSC_COMM_WORLD,"\nNumber of linear iterations: %d",lits);
					PetscPrintf(PETSC_COMM_WORLD,"\nTotal inner-products : %d", total_ip_dots);
					// PetscPrintf(PETSC_COMM_WORLD,"\nOperator applications: %d",noperator);
					// PetscPrintf(PETSC_COMM_WORLD,"\nJacobi[%D]-Preconditioned-Operator applications: %d", prec_type, prec_calls);
					PetscPrintf(PETSC_COMM_WORLD,"\n[%D]\t:Eigenvalue from SNES solve = %16.14f\n\n", iev, evals_r[nconv]);
					total_ip_dots = 0;

					which_ortho.push_back(PETSC_TRUE) ;
					nconv++ ;
				}
				// tot_operator_count += noperator ;
				// tot_prec_count += prec_calls ;
				
				// VecView(eigsol.vec(), 	PETSC_VIEWER_STDOUT_WORLD ) ;
				// std::cin.get() ;
				
			}
			
			ierr = IPDestroy(&ip_ortho) ;CHKERRQ(ierr);
			
			perf_log.pop("hybrid_newton_solve");
			
			// Optional: Get some information from the solver and display it
			PetscPrintf(PETSC_COMM_WORLD,"\n\n Solution method: Newton solve\n");
			PetscPrintf(PETSC_COMM_WORLD,"\nNumber of requested eigenvalues: %d", nev_req);
			PetscPrintf(PETSC_COMM_WORLD,"\nNumber of converged eigenvalues: %d", nconv);
			// PetscPrintf(PETSC_COMM_WORLD,"\nTotal Operator applications: %d",tot_operator_count);
			// PetscPrintf(PETSC_COMM_WORLD,"\nTotal Jacobi[%D]-Preconditioned-Operator applications: %d",prec_type, tot_prec_count);
			
			eigres.clear() ;
			eigsol.clear() ;
			
			ierr = SNESDestroy(&snes) ;CHKERRQ(ierr);
			ierr = MatDestroy(&Jac);CHKERRQ(ierr);
			ierr = MatDestroy(&Prec);CHKERRQ(ierr);
			
		}

		{
			finish = MPI_Wtime();
			run_time = finish-start;
			PetscPrintf(PETSC_COMM_WORLD, "\n\nPrimary calculation run time : %G (seconds)\n\n", run_time) ;
		}
		
		perf_log.push("compute_eigen_modes");
		
		if (nconv>0)
		{
			// Display eigenvalues and relative errors
			ierr = PetscPrintf(PETSC_COMM_WORLD,
			 "           k          ||Ax-kx||/||kx||\n"
			 "   ----------------- ------------------\n" );CHKERRQ(ierr);
			
			Vec temp; 
			ierr = VecDuplicate(evecs[0], &temp) ;CHKERRQ(ierr);
			
			PetscScalar kr, ki, error ;
			for(int iev=0; iev<nconv; iev++ )
			{
				// Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
				// ki (imaginary part)
				kr = evals_r[iev] ;		ki = evals_c[iev] ;	
				
				// ierr = VecCopy(evecs[iev], pS.vec()) ;CHKERRQ(ierr);
				
				// neutronics.Keff() = kr ;
				
				// Compute the relative error associated to each eigenpair
				double errnum, errden ;
			
				ierr = Wrapper::generalized_eigen_residual(&neutronics, evecs[iev], temp, kr) ;CHKERRQ(ierr);
				
				ierr = VecNorm(temp, NORM_2, &errnum);CHKERRQ(ierr);
				
				double res_eigval = 0 ;
				ierr = VecNorm(evecs[iev], NORM_2, &res_eigval) ;CHKERRQ(ierr);
				
				errden = res_eigval * kr ;
				
				error = errnum / errden ;

				// ierr = PetscPrintf(PETSC_COMM_WORLD,"\n errors : numerator =  %16.14f \t denominator =  %16.14f \t ratio = %16.14f \n",errnum,errden, error);
				
				{
					KSP& loss_solver = neutronics.loss_matrix_solver() ;
					
					ierr = KSPSolve(loss_solver, temp, pS.vec()) ;CHKERRQ(ierr);
					
					// ierr = Wrapper::precondition_loss_matrix(&neutronics, temp.vec(), pS.vec()) ;
					
					errnum = pS.l2_norm() ;
					
					error = errnum / errden ;

					// ierr = PetscPrintf(PETSC_COMM_WORLD,"\n errors : numerator =  %16.14f \t denominator =  %16.14f \t ratio = %16.14f \n",errnum,errden, error);
				}
				
				if (fabs(ki) > 1e-10 )	// complex eigenvalue
					ierr = PetscPrintf(PETSC_COMM_WORLD,"   %16.14f + i * %16.14f     %12.10g\n",kr,ki,error);
				else
					ierr = PetscPrintf(PETSC_COMM_WORLD,"   %16.14f       %12.10g\n",kr,error);
				CHKERRQ(ierr);
				
			}

			ierr = VecDestroy(&temp); CHKERRQ(ierr);
			
			ierr = PetscPrintf(PETSC_COMM_WORLD,"\n" );CHKERRQ(ierr);				
		}
    
		perf_log.pop("compute_eigen_modes");
		
		ierr = PetscPrintf(PETSC_COMM_WORLD,"\nMy exact keff = %16.14f\n", problem.Keff());CHKERRQ(ierr);
		
		Wrapper::clear_wrapper() ;
		
		if (iosoln)
		{
			// ierr = Wrapper::petsc_copy_vector(ctx->x, neutronics.get_solution()->vec(), false);CHKERRQ(ierr);
		
			std::stringstream out;
			perf_log.push("print_output");
			
			for (int iev=0; iev<nev_req; iev++)
			{
				// normalzie all eigenvectors
				{
					PetscScalar maxv = 0., minv = 0. ;	int ind = -1 ;
					
					ierr = VecNorm(evecs[iev], NORM_INFINITY, &maxv) ;CHKERRQ(ierr);
					ierr = VecScale(evecs[iev], 1.0/maxv) ;CHKERRQ(ierr);
					ierr = VecMax(evecs[iev], &ind, &maxv) ;CHKERRQ(ierr);
					
					if (fabs(1.0-maxv) > 1e-10)	// possibly negative maxima
					{
						ierr = VecScale(evecs[iev], -1.0) ;CHKERRQ(ierr);
					}
					
					ierr = VecMax(evecs[iev], &ind, &maxv) ;CHKERRQ(ierr);
					ierr = VecMin(evecs[iev], &ind, &minv) ;CHKERRQ(ierr);
					PetscPrintf(PETSC_COMM_WORLD, "\nMaximum value = %G and Minimum value = %G\n", maxv, minv) ;
				}
				
				ierr = VecCopy(evecs[iev], pS.vec()) ;CHKERRQ(ierr);
				neutronics.update() ;
				
				out << "problem_" << PROBLEMTYPE << "_mode_" << iev << ".vtu"  ;
				
				std::string fprefix = out.str() ;
				
				// After solving the system write the solution
				//GMVIO (neutronics.get_mesh()).write_equation_systems ("results.gmv", neutronics);
				// VTKIO (neutronics.get_mesh()).write_equation_systems (out.str(), neutronics);
				VTKIO_NOVTK (neutronics.get_mesh()).write_equation_systems (out.str(), neutronics);
				
				out.str("") ;
				
				ierr = VecDestroy(&evecs[iev]) ;CHKERRQ(ierr);
			}
			perf_log.pop("print_output");
		}
		else
		{
			for (int iev=0; iev<nconv; iev++)
			{
				ierr = VecDestroy(&evecs[iev]) ;CHKERRQ(ierr);
			}
		}
		
		evecs.clear () ;
		evals_r.clear () ;
		
		ierr = MatDestroy(&eigMatrix) ;CHKERRQ(ierr);
		
		ierr = EPSDestroy(&eps);CHKERRQ(ierr);

		Gain = PETSC_NULL ;
		Loss = PETSC_NULL ;
				
		neutronics.clear() ;
	}
	
	// All done.  
	return 0 ;	
}

/*
PetscErrorCode orthonormalize_vector(Vec soln, Vec work)
{
	PetscErrorCode ierr = 0 ;
	PetscScalar *ortho_h, *swork ;
	PetscReal ortho_norm ;
	PetscBool lindep ;
	int ndots ;
	
	// let us orthogonalize the input vector against all the previously converged vectors
#if SLEPC_VERSION_LESS_THAN(3,0,0)
	ierr = IPOrthogonalize(ip_ortho, evecs.size(), &which_ortho[0], &evecs[0], soln, 
					ortho_h, &ortho_norm, &lindep, work) ;CHKERRQ(ierr);
#else
	ierr = IPOrthogonalize(ip_ortho, evecs.size(), &which_ortho[0], &evecs[0], soln, 
					ortho_h, &ortho_norm, &lindep, work, swork) ;CHKERRQ(ierr);
#endif

	// ierr = IPView(ip_ortho, PETSC_VIEWER_STDOUT_WORLD) ;CHKERRQ(ierr);
	ierr = IPGetOperationCounters(ip_ortho, &ndots) ;CHKERRQ(ierr);
	// if (lindep == PETSC_FALSE)
		// PetscPrintf(PETSC_COMM_WORLD,"IPOrthogonalize refinement did not improve the quality of orthogonalization.\n");
	// PetscPrintf(PETSC_COMM_WORLD,"Number of inner-products during orthogonalization: %d",ndots);

	// ierr = VecScale(soln, 1.0/ortho_norm) ;
	
	total_ip_dots += ndots ;
	return ierr ;
}
*/
